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Abstract We show how the recently developed theory of geodesic transport bar- 
riers for fluid flows can be used to uncover key invariant manifolds in externally 
forced, one-degree-of-freedom mechanical systems. Speciflcally, invariant sets in 
such systems turn out to be shadowed by least-stretching geodesies of the Cauchy- 
Green strain tensor computed from the flow map of the forced mechanical system. 
This approach enables the flnite-time visualization of generalized stable and un- 
stable manifolds, attractors and generalized KAM curves under arbitrary forcing, 
when Poincare maps are not available. We illustrate these results by detailed vi- 
sualizations of the key flnite-time invariant sets of conservatively and dissipatively 
forced Duffing oscillators. 



1 Introduction 



A number of numerical and analytical techniques are available to analyze exter- 
nally forced nonlinear mechanical systems. Indeed, perturbation methods, Lya- 
punov exponents, Poincare maps, phase space embeddings and other tools have 
been become broadly used in mechanics [HH]. Still, most of these techniques, are 
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only applicable to nonlinear systems subject to autonomous (time-independent), 
time-periodic, or time-quasiperiodic forcing. 

These recurrent types of forcing allow for the analysis of asymptotic features 
based on a finite-time sample of the underlying fiow map— the mapping that takes 
initial conditions to their later states. Indeed, to understand the phase space dy- 
namics of an autonomous system, knowing the flow map over an arbitrary short 
(but finite) time interval is enough, as all trends can be reproduced by the re- 
peated applications of this short-time map. Similarly, the period map of a time- 
periodic system (or a one-parameter family of fiow maps for a time-quasiperiodic 
system) renders asymptotic conclusions about recurrent features, such as periodic 
and quasiperiodic orbits, their stable and unstable manifolds, attractors, etc. 

By contrast, the identification of key features in the response of a nonlinear 
system under time-aperiodic forcing has remained an open problem. Mathemati- 
cally, the lack of precise temporal recurrence in such systems prevents the use of a 
compact extended phase space on which the forced system would be autonomous. 
This lack of compactness, in turn, renders most techniques of nonlinear dynamics 
inapplicable. Even more importantly, a finite-time understanding of the fiow map 
can no longer be used to gain a full understanding of a (potentially ever-changing) 
non-autonomous system. 

Why would one want to develop an understanding of mechanical systems under 
aperiodic, finite-time forcing conditions? The most important reason is that most 
realistic forms of forcing will take time to build up, and hence will be transient in 
nature, at least initially. Even if the forcing is time-independent, the finite-time 
transient response of a mechanical system is often crucial to its design, as the 
largest stresses and strains invariably occur during this period. 

Similar challenges arise in fiuid dynamics, where temporally aperiodic unsteady 
fiows are the rule rather than the exception. Observational or numerical data for 
such fluid flows is only available for a limited time interval, and some key features of 
the flow may only be present for an even shorter time. For instance, the conditions 
creating a hurricane in the atmosphere are transient, rather than periodic, in 
nature, and the hurricane itself will generally only exist for less than two weeks 
[3]. As a result, available asymptotic methods are clearly inapplicable to its study, 
even though there is great interest in uncovering its internal structure and overall 
dynamics. 

In response to these challenges in fluid dynamics, a number of diagnostic tools 
have been developed [35]. Only very recently, however, has a rigorous mathemat- 
ical theory emerged for dynamical structures in finite-time aperiodic flow data [6] . 
This theory flnds that finite-time invariant structures in a dynamical system are 
governed by intrinsic, metric properties of the flnite-time flow map. Speciflcally, 
in two-dimensional unsteady flows, structures acting as transport barriers can be 
uncovered with the help of geodesies of the Cauchy-Green strain tensor used in 
continuum mechanics [7]- This approach generalizes and extends earlier work on 
hyperbolic Lagrangian Coherent Structures (LCS), that are locally most repelling 
or attracting material lines in the flow [8ll9l ll0l[TT| . 

In this paper, we review the geodesic transport theory developed in [6] in 
the context of one-degree-of-freedom, aperiodically forced mechanical systems. We 
then show how this theory uncovers key invariant sets under both conservative and 
dissipative forcing in cases where classic techniques, such as Poincare maps, are not 
available. Remarkably, these flnite-time invariant sets can be explicitly identified 



Title Suppressed Due to Excessive Length 



3 



as parametrized curves, as opposed to plots requiring post-processing or feature 
extraction. 

The organization of this paper is as foUows. Section j|2]is divided into two sub- 
sections: Section |2.1| provides the necessary background for the geodesic theory of 
transport barriers developed in [6] . In section ^2.2| we describe a numerical imple- 
mentation of this theory that detects finite-time invariant sets as transport bar- 
rier. Section !|3] presents results from the application of this numerical algorithm 
to one degree-of-freedom mechanical systems. First, as a proof of concept, |3.1| 
considers conservative and dissipative time-periodic Duffing oscillators, compar- 
ing their geodesically extracted invariant sets with those obtained form Poincare 
maps. Next, section j ]3.2| deals with invariant sets in aperiodically forced Duffing 
oscillators, for which Poincare maps or other rigorous extraction methods are not 
available. We conclude the paper with a summary and outlook. 



2 Set-up 

The key invariant sets of autonomous and time-periodic dynamical systems-such 
as fixed points, periodic and quasiperiodic motions, their stable and unstable man- 
ifolds, and attractors-are typically distinguished by their asymptotic properties. 
In contrast, invariant sets in finite-time, aperiodic dynamical systems solely dis- 
tinguish themselves by their observed impact on trajectory patterns over the finite 
time interval of their definition. This observed impact is a pronounced lack of tra- 
jectory exchange (or transport) across the invariant set, which remains coherent 
in time, i.e., only undergoes minor deformation. Well-understood, classic exam- 
ples of such transport barriers include local stable manifolds of saddles, parallel 
shear jets, and KAM tori of time-periodic conservative systems. Until recently, 
a common dynamical feature of these barriers has not been identified, hindering 
the unified detection of transport barriers in general non-autonomous dynamical 
systems. 

As noted recently in [6J, however, a common feature of all canonical transport 
barriers in two dimensions is that they stretch less under the flow than neighboring 
curves of initial conditions do. This observation leads to a nonstandard calculus 
of variations problem with unknown endpoints and a singular Lagrangian. Below 
we recall the solution of this problem from "6^ , with a notation and terminology 
adapted to one-degree-of-freedom mechanical oscillators. 

A one-degree-of-freedom forced nonlinear oscillator can generally be written as 
a two-dimensional dynamical system 

± = v{x,t), xeUcR^, t£[to,ti], (1) 

with U denoting an open set in the state space, where the vector x labels tuples of 
positions and velocities. The vector v{x,t), assumed twice continuously differen- 
tiable, contains the velocity and acceleration of the system at state x and at time 
t. 

Let x(ti;to,xo) denote the final state of system ^ at time ti, given its state 
xq at an initial time to. The fiow map associated with ^ over this time interval 
is defined as 

:xo^x(ti;to,a^o), (2) 
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which maps initial states to final states at ti. The Cauchy-Green (CG) strain 
tensor associated with the flow map ([2| is defined as 

{xo) = [DF}^ {xo)VdfI^ (xo), (3) 

where DF^^ denotes the gradient of the fiow map and the symbol T refers to 
matrix transposition. 

Note that the CG tensor is symmetric and positive definite. As a result, it has 
two positive eigenvalues < Ai < A2 and an orthonormal eigenbasis {^i,i^2}- We 
fix this eigenbasis so that 

Cto(a;o)?i(2;o) = Ai(a;o)Ci(2;o), |Ci(a;o)| = 1, i€{l,2}, 

U^o) = o^iixo), n=(^^^^^y (4) 

We suppress the dependence of and on to and ti for notational simplicity. 



2.1 Geodesic transport barriers in phase space 

A material line 74 = F^^ (74^, ) is an evolving curve of initial conditions jto under 
the fiow map F^^. As shown in [6], for such a material line to be a locally least- 
stretching curve over [to,ii], it must be a hyperbolic, a parabolic or an elliptic line 
(see figure [TJ . 

The initial position 74^ of a hyperbolic material line is tangent to the vector 
field ^1 at all its points. Such material lines are compressed by the fiow by locally 
the largest rate, while repelling all nearby material lines at an exponential-in-time 
rate. The classic example of a hyperbolic material lines is the unstable manifold 
of a saddle- type fixed point. 

A parabolic material line is an open material curve whose initial position 74^ is 
tangent to one of the directions of locally largest shear. At each point of the phase 
space, the two directions of locally largest shear are given by 



as derived in [6]. Parabolic material lines still repel most nearby material lines 
(except for those parallel to them), but only at a rate that is linear in time. Classic 
examples of parabolic material lines in fiuid mechanics are the parallel trajectories 
of a steady shear fiow. 

Finally, an elliptic material line is a closed curve whose initial position 74^ is 
tangent to one of the two directions of locally largest shear given in ^. As a 
result, elliptic lines also repel nearby, nonparallel material lines at a linear rate, 
but they also enclose a connected region. Classic examples of elliptic material lines 
are closed trajectories of a steady, circular shear flow, such as a vortex. 

Initial positions of hyperbolic material lines are, by deflnition, stramlines, i.e., 
trajectories of the autonomous differential equation 



r' = ^i(r), ref/cR^ 



(6) 
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(b) A parabolic barrier (red curve) is an open curve that has the locally largest rate of 
Lagrangian shear along its tangent. 



t^to t = ti 




(c) An elliptic barrier (red curve) is a closed curve with the same dynamical property as 
a paraolic barrier. 

Fig. 1: The three types of transport barriers in two-dimensional flows. 



where r : [0,1] i-^ U is the parametrization of the strainhne by arc-length. A 
hyperbolic barrier is then a strainhne that is locahy the closest to least-stretching 
geodesies of the CG tensor, with the latter viewed as a metric tensor on the domain 
U of the phase space. The pointwise closeness of strainlines to least-stretching 
geodesies can be computed in terms of the invariants of the CG strain tensor. 
Specifically, the distance (difference of tangents plus difference of curvatures) 
of a strainhne from the least-stretching geodesic of Cl^ through a point xq is given 
by the geodesic strain deviation 
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,,.(^^)^KvW|Lp^^ (7) 

with Ki(a;o) denoting the curvature of the strainhne through xo [6]. A hyperbohc 
barrier is a compact strainhne segment on which is pointwise below a smaU 
threshold value, and whose averaged value is locally minimal relative to all 
neighboring strainlines. 

Similarly, initial positions of parabolic and elliptic material lines are, by defi- 
nition, shearlines, i.e., trajectories of the autonomous differential equation 

r' = r?±(r), r e U C . (8) 

A parabolic barrier is an open shearline that is close to least-stretching geodesies 
of the CG tensor. The pointwise C^-closeness of shearlines to least-stretching 
geodesies is given by the geodesic shear deviation 



(VA2,Ci) (VA2,C2>(^/^+A^'-^A2') 



A^i [Vm"" + (1 - Ai) VTTM 



AlVTTAa VTTAz' 



+ (9) 



with k,2{xq) denoting the curvature of the ^2 vector field at the point xq 6 . The 
geodesic shear deviation should pointwise be below a small threshold level for an 
open shearline to qualify as a parabolic barrier. Similarly, a closed shearline is 
an elliptic barrier if its pointwise geodesic shear deviation is smaller than small 
threshold level. 

For the purposes of the present discussion, we call a mechanical system of the 
form ([1]) conservative if it has vanishing divergence, i.e., V ■ v{x,t) = 0, with V 
referring to differentiation with respect to x. This property implies that flow map 
of ([1]) conserves phase-space area for all times [13] . 

While a typical material line in such a conservative system will still stretch and 
deform significantly over time, the length of a shearline will always be preserved 
under the area- preserving flow map (cf. [6]). An elliptic barrier in a conser- 
vative system will, therefore, have the same enclosed area and arclength at the 
initial time Iq and at the final time t\. These two conservation properties imply 
that an elliptic barrier in a non-autonomous conservative system may only undergo 
translation, rotation and some slight deformation, but will otherwise preserve its 
overall shape. As a result, the interior of an elliptic barrier will not mix with the 
rest of the phase-space, making elliptic barriers the ideal generalized KAM curves 
in aperiodically forced conservative mechanical systems. 



2.2 Computation of invariant sets as transport barriers 

In this section, we describe numerical algorithms for the extraction of hyperbolic 
and elliptic barriers in a one-degree-of-freedom mechanical system with general 
time dependence. Parabolic barriers can in principle also exist in mechanical sys- 
tems, but they do not arise in the simple examples we study below. In contrast, 
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parabolic barriers are more common in geophysical fluid mechanics where they 
typically represent unsteady shear jets. 

Our numerical algorithms require a careful computation of the CG tensor. In 
most mechanical systems, trajectories separate rapidly, resulting in an exponential 
growth in the entries of the CG tensor. This growth necessitates the use of a well- 
resolved grid, as well as the deployment of high-end integrators in solving for the 
trajectories of ([l]) starting form this grid. Further computational challenges arise 
from the handling of the unavoidable orientational discontinuities and isolated 
singularities of the eigenvector fields and ^2- The reader is referred to Farazmand 
& Haller [10] for a detailed treatment of these computational aspects. 

As a zeroth step, we fix a sufficiently dense grid Qo of initial conditions in the 
hase-space U, then advect the grid points from time to to time ti under system 



jj). This gives a numerical representation of the flow map F^^ over the grid Qo- 
The CG tensor field Cj^ is then obtained by definition tih from Ff^ . In computing 



the gradient F>F*^ , we use careful finite differencing over an auxiliary grid, as 
described in [10'. 

Since, at each point xo & Qo, the tensor Cl^{xo) is a two- by-two matrix, com- 
puting its eigenvalues {Ai,A2} and eigenvectors {^i,.f2} is straightforward. With 
the CG eigenvalues and eigenvectors at hand, we locate the hyperbolic barriers 
using the following algorithm. 

Algorithm 1 (Locating hyperbolic barriers) 

1. Fix a small positive parameter e^^ as the admissible upper bound for the point- 
wise geodesic strain deviation of hyperbolic transport barriers. 

2. Calculate strainlines by solving the ODE ([6| numerically, with linear interpo- 
lation of the strain vector field between grid points. Truncate strainlines to 
compact segments whose pointwise geodesic strain deviation is below e^^ 

3. Locate hyperbolic barriers as strainline segments jto with locally minimal rel- 
ative stretching, i.e., strainline segments that locally minimize the function 

(10) 

Here l{'yto) and l{'yti) denote the length of the strainline 7t(, and the length of 
its advected image jti , respectively. 



Computing the relative stretching ( |10[ ) of a strainline 'yt^, in principle, requires 
advecting the strainline to time ti. However, as shown in [6 , the length of the 
advected image satisfies i(7ti) = f where the integration is carried out 

along the strainline 7tp . This renders the strainline advection unnecessary. 

Numerical experiments have shown that a direct computation of is usually 
less accurate than that of ^2 due to the attracting nature of strongest eigenvector 
of the CG tensor [W. For this reason, computing ^1 as an orthogonal rotation of 
^2 is preferable. Moreover, it has been shown ,12_ that strainlines can be computed 
more accurately as advected images of stretchlines, i.e. curves that are everywhere 
tangent to the second eigenvector of the backward-time CG tensor Cj° . In the 
present paper, this approach is taken for computing the strainlines. 

Computing elliptic barriers amounts to finding limit cycles of the ODE ([8|. 
To this end, we follow the approach used in [6lll2| by first identifying candidate 
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Poincare Section 




Fig. 2: Locating closed shearlines using a Poincare section of the shear vector field. 
Closed shearlines pass through the fixed points of the corresponding Poincare map. 



regions for shear limit cycles visually, then calculating the Poincare map on a one- 
dimensional section transverse to the flow within the candidate region (see figure 
[2| . Hyperbolic fixed points of this map can be located by iteration, marking limit 
cycles of the shear vector field (see [12] for more detail). 

This process is used in the following algorithm to locate elliptic barriers. 

Algorithm 2 (Locating elliptic barriers) 

1. Fix a small positive parameter Srj^ as the admissible upper bound for the 
average geodesic shear deviation of elliptic transport barriers. 

2. Visually locate the regions where closed shearlines may exist. Construct a suf- 
ficiently dense Poincare map, as discussed above. Locate the fixed points of the 
Poincare map by iteration. 

3. Compute the full closed shearlines emanating from the fixed points of the 
Poincare map. 

4. Locate elliptic barriers as closed shearlines whose average geodesic deviation 
(dg* ) satisfies (dg* ) < . 

In the next section, we use the above algorithms for locating invariant sets in 
simple forced and damped nonlinear oscillators. 



3 Results 

We demonstrate the implementation of the geodesic theory of transport barriers 
on four Duffing-type oscillators. As a proof of concept, in the first two examples 
(section S3.ll, we consider periodically forced Duffing oscillators for which we can 
explicitly verify our results using an appropriately defined Poincare map. 
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The next two examples deal with aperiodically forced Duffing oscillators (section 



!3.2). In these examples, despite the absence of a Poincare map, we still obtain 
the key invariant sets as hyperbolic and elliptic barriers. 

To implement algorithms 1 and 2 in the forthcoming examples, the CG tensor 
is computed over a uniform grid Qo of 1000 x 1000 points. A fourth order Runge- 
Kutta method with variable step-size (ODE45 in MATLAB) is used to solve the 
first-order ODEs ([T]), ([6| and ^ numerically. The absolute and relative tolerances 
of the ODE solver are set equal to 10^"' and 10^®, respectively. Off the grid points, 
the strain and shear vector fields are obtained by bilinear interpolation. 

In each case, the Poincare map of algorithm [2] is approximated by 500 points 
along the Poincare section. The zeros of the map are located by a standard secant 
method. 



3.1 Proof of concept: Periodically forced Duffing oscillator 
Case 1: Pure periodic forcing, no damping 
Consider the periodically forced Duffing oscillator 

X2 = xi — Xi + ecos(t). 

For e = 0, the system is integrable with one hyperbolic fixed point at (0, 0), and 
two elliptic fixed points (1,0) and (—1,0), respectively. As is well known, there are 
two homoclinic orbits connected to the hyperbolic fixed point, each enclosing an 
elliptic fixed point, which is in turn surrounded by periodic orbits. These periodic 
orbits appear as closed invariant curves for the Poincare map P := Fq^ . The fixed 
points of the fiow are also fixed points of P. 

For < e <C 1, the Kolmogorov-Arnold-Moser (KAM) theory |13] guarantees 
the survival of most closed invariant sets for P. Figure |3] shows these surviving 
invariant sets (KAM curves) of P obtained for e = 0.08. For the KAM curves 
to appear continuous-looking, nearly 500 iterations of P were needed, requiring 
the advection of initial conditions up to time t = IOOOtt. The stochastic region 
surrounding the KAM curves is due to chaotic dynamics arising from the transverse 
intersections of the stable and unstable manifold of the perturbed hyperbolic fixed 
point of P. 

The surviving KAM curves are well-known, classic examples of transport bar- 
riers. We would like to capture as many of them as possible as elliptic barriers 
using the geodesic transport theory described in previous sections. Note that not 
all KAM curves are expected to prevail as locally least-stretching curves for a 
given choice of the observational time interval [to,ti]; some of these curves may 
take longer to prevail due to their shape and shearing properties. 



We use the elliptic barrier extraction algorithm of section 2.2 with e,,^ = 0.7. 
Figure [4] shows the resulting shearlines in the KAM regions, with the closed ones 
marked by red. Note that these shearlines were obtained from the CG tensor com- 
puted over the time interval [0,87r], spanning just four iterations of the Poincare 
map. Despite this low number of iterations, the highlighted elliptic barriers are 
practically indistinguishable form the KAM curves obtained from five hundred 
iterations. 
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Fig. 3: Five hundred iterations of the Poincare map for the periodicaUy forced 
Duffing osciUator. Two elhptic regions of the phase-space fiUed by KAM tori are 
shown. 




Fig. 4: Shearhnes (black) of the periodically forced Duffing oscillator computed 
at to = 0, with integration time T = Stt. The extracted elliptic barriers with 
(c'g*) ^ £ri± =0.7 are shown in red. 



Figure [5] shows the convergence of an elliptic barrier to a KAM curve as the 
integration time T = ti — to increases. Note how the average geodesic deviation 
(dg*) decreases with increasing T, indicating decreasing deviation from nearby 
Cauchy-Green geodesies. 

Remarkably, constructing these elliptic barriers requires significantly shorter 
integration time (only four forcing periods) in comparison to visualization through 
the Poincare map, which required 500 forcing periods to reveal KAM curves as 
continuous objects. Clearly, the overall computational cost for constructing elliptic 
barriers still comes out to be higher, since the CG tensor needs to be constructed 
on a relatively dense grid Go, a-s discussed in section §2.2[ This high computational 



cost will be justified, however, in the case of aperiodic forcing (section [ 3.2 1, where 
no Poincare map is available. 

In the context of one-degree-of-freedom mechanical systems, the outermost el- 
liptic barrier marks the boundary between regions of chaotic dynamics and regions 
of oscillations that are regular on a macroscopic scale. To demonstrate this sharp 
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_Q4l , , , , , _Q4l , , , , , _g4l , , , , , 

O.a 0.9 1 1.1 1.2 1.3 O.B 0.9 1 1.1 1.2 1.3 0.8 0.9 1 1.1 1.2 1. 



(a) T = Svr, (4±> = 0.5991. (b) T = lOvr, (4±> = 0.3889. (c) T = 127r, (4±> = 0.2355. 

Fig. 5: Convergence of an elliptic barrier (red) to a KAM curve (black) as the 
integration time T = ti ~ to increases. The gradually decreasing average geodesic 
deviation (dg* ) confirms the convergence to Cauchy-Green geodesies that closely 
shadow the underlying KAM torus. 



dividing property of elliptic barriers, we show the evolution of system (12 1 from 
three initial states, two of which are inside the elliptic region and one of which is 
outside (figure |6^). The system exhibits rapid changes in its state when started 
from outside the elliptic region. In contrast, more regular behavior is observed for 
trajectories starting inside the elliptic region. This behavior is further depicted in 
figure [Tj which shows the evolution of the xi-coordinate of the trajectories as a 
function of time. 

Case 2: Periodic forcing and damping 

Consider now the damped-forced Dufhng oscillator 

il = X2, 

X2 = xi — xi — 5x2 + ecos(t), (11) 

with S = 0.15 and e = 0.3. This system is known to have a chaotic attractor 
that appears as an invariant set of the the Poincare map P = Fg^ (see, e.g., [1]). 
Here, we show that the attractor can be very closely approximated by hyperbolic 
barriers computed via algorithm [l] 

Figure shows strainlines computed backward in time with to = and in- 
tegration time T = ti — to = — Svr. The strainline with globally minimal relative 



stretching ( 10 1 is shown in figure|8p. Black dots mark the points where the geodesic 
deviation d^' exceeds the admissible upper bound e^^ = 10~^. At its tail (covered 
by black dots), the strainline persistently deviates from CG geodesies, and hence 
should be truncated. The resulting hyperbolic barrier, as a finite-time approxima- 
tion to the chaotic attractor, is shown in figure [8]:. 

The approximate location of the attractor can also be revealed by applying 
the Poincare map to a few initial conditions (tracers) released from the basin of 
attraction. For long enough advection time, the initial conditions converge to the 
attractor highlighting its position (see figure[9|i and[9j3). In figure|9]:, the hyperbolic 
barrier is superimposed on the advected tracers showing close agreement between 
the two. Figure^ shows the tracers advected for a longer time (T = AQtv) together 
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Fig. 6: (a) The outermost elliptic barrier (black curve) and three initial conditions: 
Two inside the elliptic barrier (blue and green) and one outside the elliptic barrier 
(red), (b) The corresponding trajectories are shown in the extended phase space of 
{xi, X2,t). The closed black curves mark the elliptic barrier at to = and ti = IGtt. 



1.5 
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-2 -1.5 -1 -0.5 0.5 1 1.5 2 
(a) Strainlines co mput ed for the damped-forced 
Duffing oscillator l |ll| at time to = 0, with the 
integration time T = — Stt. 





-1.5 -1 -0.5 0.5 1 1.5 

(b) The strainline (red) with globally minimum 



highlighted as black dots. 



-1.5 -1 -0.5 0.5 1 1.5 

(c) Final approximation of the chaotic attractor 
by a single, continuous strainline with minimal 
geodesic deviation. 



Fig. 8: Construction of the attractor of the damped-forced Duffing osciUator as a 
hyperbolic transport barrier. 



with the hyperbohc barrier; the two virtuaUy coincide. Note that the hyperbohc 
barrier is a smooth, parametrized curve (computed as a trajectory of (|6|), while 
the tracers form a set of scattered points. 



3.2 The aperiodicaUy forced Duffing osciUator 



In the next two examples, we study aperiodicaUy forced Duffing oscillators. In the 
presence of aperiodic forcing, the Poincare map P is no longer defined as the sys- 
tem lacks any recurrent behavior. However, KAM-type curves (i.e., closed curves, 
resisting significant deformation) and generalized stable and unstable manifolds 
(i.e., most repelling and attracting material lines) exist in the phase-space and 
determine the overall dynamics of the system. 
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(a) 



(b) 





(c) 



(d) 



Fig. 9: (a) Attractor of system (11) obtained from four iterates of the Poincare 
map. (b) Attractor obtained from 20 iterates of the Poincare map. (c) Attractor 
computed as a hyperboHc barrier (red), compared with the Poincare map (blue) 
computed for the same integration time (four iterates), (d) Comparison of attractor 
computed as a hyperbolic barrier (red) with the one obtained from 20 iteration of 
the Poincare map (blue). The integration time for locating the hyperbolic barrier 
is T = ti - to = -Stt. 



Case 1: Purely aperiodic forcing, no damping 
Consider the DufHng oscillator 

il = X2, 

±2=xi^xl + f{t), (12) 

where f{t) is an aperiodic forcing function obtained from a chaotic one-dimensional 
map (see figure 10 1. 

While, KAM theory is no longer applicable, one may still expect KAM-type 
barriers to survive for small forcing amplitudes. Such barriers would no longer be 
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Fig. 10: Chaotic forcing func- 



tion f{t) for equation ( 12 ) 



repeating themselves periodically in the extended phase space. Instead, a gener- 
alized KAM barrier is expected to be an invariant cylinder, with cross sections 
showing only minor deformation. The existence of such structures can, however, 
be no longer studied via Poincare maps. 

Figure [TT] confirms that generalized KAM- type curves, obtained as elliptic bar- 
riers, do exist in this problem. These barriers are computed over the time interval 
[0, 47r] (i.e. to = and ti = to+T = Aiv). As discussed in section { 2.1 the arclength 



of an elliptic barrier at the initial time to is equal to the arclength of its adverted 
image under the fiow map F^^ at the final time ti. This arclength preservation is 
illustrated numerically in figure [121 which shows the relative stretching. 



Sl{t) 



^(70) 



(13) 



of the time-t image 74 of an elliptic barrier 70, with £ referring to the arclength of 
the curve. Ideally, the relative stretching of each elliptic barrier should be zero at 
time ti = An, i.e. S£{4:Tv) = 0. Instead, we find that the relative stretching SE^A-k) 
of the computed elliptic barriers is at most 1.5%. This deviation from zero arises 



from numerical errors in the computation of the CG strain tensor C^' 
turn causes small inaccuracies in the computation of closed shearlines. 



which in 



As noted earlier, the small relative stretching and the conservation of enclosed 
area for an elliptic barrier in incompressible flow only allows for small deformations 
when the barrier is adverted in time. This is illustrated in figure [l3j which shows 
the blue elliptic barrier of figure [TT[3 in the extended phase-space. Each constant- 
time slice of the figure is the adverted image of the barrier. 
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Fig. 11: Closed shearlines for equation (121 computed in two elliptic regions. The 



figure shows the shearlines at time to = 0. The integration time is T = Att. 
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The 



Fig. 12: The relative stretching S£{t) x 100 of closed shearlines of figure 
colors correspond to those of figure [TT] By their arc-length preservation property, 
the adverted elliptic barriers must theoretically have the same arclength at times 
to = and ti = A-k. The numerical error in arclength conservation is small overall, 
but more noticeable for oscillations with large amplitudes (green and red curves 
of the right panel). 




Fig. 13: Generalized KAM-type cylinder in the extended phase space of the aperi- 
odically forced Duffing undamped oscillator. The cylinder is obtained by advection 
of the closed shearline shown in blue in figure lib). 
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Finally, we point out that the stability of the trajectories inside elliptic barriers 
show a similar trend as in the case of the periodically forced Duffing equation 
(figures [6] and [7]). Namely, perturbations inside the elliptic regions remain small 
while they grow significantly inside the hyperbolic regions. 

Case 2: Aperiodic forcing with damping 

In this final example, we consider the aperiodically forced, damped Duffing oscil- 
lator 

±1 = X2, 

±2=xi~xl~5x2 + f{t), (14) 

with damping coefficient 5 = 0.15. The forcing function f{t) is similar to that of 
Case I above, but with an amplitude twice as large. As a result, none of the elliptic 
barriers survive even in the absence of damping. 

Again, because of the aperiodic forcing, the behavior of this system is a priori 
unknown and cannot be explored using Poincare maps. In order to investigate the 



existence of an attractor, strainlines (figure 14 1) are computed from the backward- 



time CG strain tensor Cf ' with tn = 30 and ti = 10. The strainline with minimum 



relative stretching ( 10 ) is then extracted. The part of this strainline satisfying 
is considered as the most influential hyperbohc barrier (figure h4b). The 



admissible upper bound e^^ for the geodesic deviation is fixed as 10^^. 

In order to confirm the existence of the extracted attractor, we advect tracer 
particles in forward time, first from time ti = 10 to time to = 30, then from ti = 
to time to = 30. Because of the fast-varying dynamics and weak dissipation, a 
relatively long advection time is required for the tracers to converge to the attrac- 



tor. Figure 15 shows the evolution of tracers over [ti,to]- Note that the attractor 
inferred from the tracers is less well pronounced than the hyperbolic barrier ex- 
tracted over the same length of time. This shows a clear advantage for geodesic 
transport theory over simple numerical experiments with tracer advection. For a 
longer integration time from to = to t = 30, the tracers eventually converge to 
the hyperbolic barrier. 

Repelling hyperbolic barriers can be computed similarly using forward-time 



computations. Figure 16 shows both hyperbolic barriers (stable and unstable man- 
ifolds) at time to = 30. The repelling barrier is computed from the CG strain tensor 
C*J with to = 30 and ti = 50. 



4 Summary and Conclusions 

We have shown how the recently developed geodesic theory of transport barriers 
in fiuid fiows can be adapted to compute finite-time invariant sets in one-degree- 
of-freedom mechanical systems with general forcing. Specifically, in the presence 
of general time dependence, temporally aperiodic stable- and unstable manifolds, 
attractors, as well as generalized KAM tori can be located as hyperbolic and ellip- 
tic barriers, respectively. The hyperbolic barriers are computed as distinguished 
strainlines, i.e. material lines along which the Lagrangian strain is locally maxi- 
mized. The elliptic barriers, on the other hand, appear as distinguished shearlines, 
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i.e. material lines along which the Lagrangian shear is locally maximized. The 
barriers are finally identified as strainlines and shearlines that are most closely ap- 
proximated by least-stretching geodesies of the metric induced by Cauchy-Green 
strain tensor. 

We have used four simple examples for illustration. First, as benchmarks, we 
considered periodically forced Duffing equations for which stable and unstable 
manifolds, attractors and KAM curves can also be obtained as invariant sets of an 
appropriately defined Poincare map. We have shown that elliptic barriers, com- 
puted as closed shearlines, coincide with the KAM curves. Also, stable and un- 
stable manifolds, as well as attractors, can be recovered as hyperbolic barriers. 
More precisely, as the integration time T = ti — to of the Cauchy-Green strain 
tensor C^J increases, the elliptic barriers in the periodically forced Duffing equa- 
tions converge to KAM curves. Similarly, the chaotic attractor of the periodically 
forced and damped Duffing equation is more and more closely delineated by a 
hyperbolic barrier computed from the backward-time Cauchy-Green strain tensor 
Cj° for increasing T = to — ti where to > ti. 

In the second set of examples, we have computed similar structures for an aperi- 
odically forced Duffing oscillator with and without damping. In this case, Poincare 
maps are no longer well-defined for the system, and hence we had to advect tracer 
particles to verify the predictions of the geodesic theory. Notably, tracer advection 
takes longer time to reveal the structures in full detail than the geodesic theory 
does. Also, tracer advection is only affective as a visualization tool if it relies on a 
small number of particles, which in turn assumes that one already roughly knows 
the location of the invariant set to be visualized. Finally, unlike scattered tracer 
points, geodesic barriers are recovered as parametrized smooth curves that provide 
a solid foundation for further analysis or highly accurate advection. 

In our examples, elliptic barriers have shown themselves as borders of subsets 
of the phase-space that barely deform over time. In fact, as illustrated in figure 
|6] outermost elliptic barriers define the boundary between chaotic and regular 
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Fig. 15: (a) Tracers advected over the time interval from ti = 10 to to = 30. 
(b) Tracers advected over a longer time interval from ti = to to = 30. (c) The 
hyperbolic barrier (red) superimposed on the tracers advected for the same time 
interval (d) Comparison of the hyperbolic barrier (red) with the tracers advected 
for the longer time interval. 



dynamics. Trajectories initiated inside elliptic barriers remain confined and ro- 
bust with respect to small perturbations. We believe that this property could be 
exploited for stabilizing mechanical systems with general time dependence. For 
instance, formulating an optimal control problem for generating elliptic behavior 
in a desired part of the phase-space is a possible approach. 

Undoubtedly, the efficient and accurate computation of invariant sets as geodesic 
transport barriers requires dedicated computational resources. Smart algorithms 
reducing the computational cost are clearly of interest. Parallel programming (both 
at CPU and GPU levels) has previously been employed for Lagrangian coherent 
structure calculations and should be useful in the present setting as well (see e.g. 
[14]). Other adaptive techniques are also available to lower the numerical cost by 
reducing the computations to regions of interest (see e.g. [151116] ). 
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Fig. 16: Attracting (blue) and repelling (red) barriers at to = 30 extracted from 
backward-time and forward-time computations, respectively. 

In principle, invariant sets in higher-degree-of-freedom mechanical systems 
could also be captured by similar techniques as locally least-stretching surfaces. 
The development of the underlying multi-dimensional theory and computational 
platform, however, is still underway. 
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